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Abstract — The design of deterministic filters can be cast 
as a problem of minimizing an associated cost function 
for an optimal control problem. Employing the min-plus 
linearity property of the dynamic programming operator 
(associated with the control problem) results in a compu- 
tationally feasible approach (while avoiding linearization 
of the system dynamics/output). This article describes the 
salient features of this approach and a specific form of 
pruning/projection, based on clustering, which serves to 
facilitate the numerical efficiency of these methods. 

I. Introduction 

Deterministic filtering approaches [1], [2], [3], [4] 
have been studied as an ahernative to stochastic methods 
of filtering. These deterministic methods are especially 
appealing in cases where the disturbance statistics are 
not known apriori. In fact, this approach has been suc- 
cessfully applied across various domains - quantitative 
finance, attitude estimation [5], etc. It is interesting to 
note that in the case of linear systems with Gaussian 
white noise disturbances/measurement noise, both these 
approaches yield the same solution (the Kalman filter 
obtained as the solution to the associated Riccati equa- 
tion): a required feature of any optimal filtering method. 
The design of deterministic filters proceeds by casting 
the filtering problem as a optimal control problem on the 
system with dynamics that are time reversed (with re- 
spect to the original system). Solving the optimal control 
problem using the dynamic programming method gives 
rise to an associated partial differential equation (the 
Hamilton-Jacobi-Bellman (HJB) equation). For systems 
with nonlinear dynamics/output equations most estima- 
tion schemes proceed via using a local linearized ap- 
proximation model. Unfortunately, issues arise in the use 
of such methods for systems with larger nonlinearities 
(c.f. [6]). A specific class of methods that have emerged 
as a promising approach to optimal controller/filter de- 
sign for nonlinear systems - the idempotent methods. 
Also termed, min/max plus methods these approaches 
exploit the fact that the dynamic programming operator 
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in the HJB equation is a linear operator in a space specif- 
ically chosen for this property [7], [8]. By constructing 
a basis expansion for the value function in this space 
(semi-field), the solution to the optimal control/filtering 
problem is rendered numerically tractable for previously 
more difficult classes of problems. These ideas were 
introduced for filtering in nonlinear systems in [8] where 
the basis used are in the (semi-convex) dual space 
(obtained via the Fenchel transform). In that work, the 
value function was transformed into this space via the 
Fenchel transform and a fixed (albeit possible countably 
infinite) set of basis functions were chosen in this dual 
space. This is in contrast to the approach herein where 
we use a set of convex functions that are not fixed 
across different time steps. More recently [9] developed 
a related approach using a different representation in 
terms of semi-convex function basis. This was termed 
curse of dimensionality free methods and have been 
utilized in areas such as quantum control [10], deception 
games [11]. Recent work [12] introduced the application 
of these methods in the context of deterministic filtering 
for nonlinear systems. Such idempotent methods have 
also been applied in other areas [13]. In this article 
we describe the design/implementation of the min-plus 
technique for deterministic filter design. Of specific 
novelty, are the approaches used to help handle the 
growth in the number of basis elements used to represent 
the value function as new measurements are made and 
state estimates are updated. 



The outline of the article is as follows. In Sec. |II] 
we introduce the problem of interest and proceed to 
describe (in Sec.HlIb the various stages of these methods: 
(i) min-plus basis expansion, (ii) recursive update of 
basis expansion in response to new information, and 
(iii) pruning of these quadratic basis terms to manage 
the computational burden of the grown in their number 
These steps are then applied to an example problem in 
Sec. HV] We conclude in Sec. [V] with an indication of 
future research directions. 



II. Problem Formulation 
Consider a system described by 

Xk+i = A{xk) + Bwh, Vk =C{xk) +Vk. 

In order to design a filter for this system we use the 
following form of the dynamics (backward dynamics 
equation) 



Xk = A{xk+i) + Bwk+i 



(1) 



and the output equation remains unchanged. Given an 
initial state estimate xq- the filter is obtained by mini- 
mizing the cost function 

1, 



JT(i,W(.)) = -||xo -xo||^ + 



fc=0 



where N, R, Q,f are the weights on the different terms 
of interest and Vk = yk — C{xk)- The nonlinear optimal 
control problem for the system in dl} corresponds to the 
following optimal cost function 

Vt{x) = vai Jt{x,wi.)). 

IU(.) 

Applying the dynamic programming approach to this 
problem leads to the following relation between the 
value function at consecutive time steps 

Vk+i{x) = mmlVkixik - l)\x{k) = x)) 

+^wo^Q„wo + -\\y-C{x)\\lY (2) 

where Vk{x{k — l)|a:;(fc) — x)) denotes the value 
function at time fc — 1 given a state x at time k. At 
T = 



Vo{x) 



\x - xa\\%o + (jP 



which can be written in the quadratic form 



A^ 
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where 



Q^i 
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(3) 
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-SjiVO, (/)? = x^N°xo + (l3" andZo = {1}. From 
[12] we recall the following result 



Theorem 2.1: Assuming that: the value function Vq 
has the quadratic form (|3); that there exist J', C such 
that 
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(4) 



and that for all M,M there exist Qg, Q" such that the 
following expansions hold 



A{xfAIA{x):^ f\ { x^ 1 ) QS(M) 

aeio 



MAix):^ /\ {x'^ 1 )QS(M) 
then, there exist /i , QZ' such that 



V,{x)^ A (^"^ 1 )Ql 
keXi 
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(5) 



(6) 



D 



The optimal state estimate at any time step k is given 



by 



X* = argmin Vt (x) . 



(7) 



It is of interest to note that x* is in fact the argmin of 
one of the quadratics in the expansion of Vk{x). Hence, 
given a form 



1 



[x^l] 
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the minimizing x* for this quadratic is given by 

X* = -[qil+qn]-^[qi2 + qIi]- (8) 



In case the states are constrained to a specific set, 
then the minimization in (|7]) must be performed in the 
permissible set of states (which would in turn change 
®). 

The above result (Thm. 12. It thus provides a recursion 
which can be applied repeatedly to determine/update 
the state estimate at every time step. In the following 
section, we describe the various stages involved in this 
approach. 

III. The stages of the min-plus approach 

At each time instant, the approach introduced above 
involves carrying out three steps: 

• Using the current value of the output to obtain a 

min-plus expansion of the output (and associated) 

functions (IH, the dynamics ^. 



• Utilizing the recursion described in Thm. 12.11 to 
obtain the new min-plus basis expansion and the 
state estimate. 

• Projection/ Pruning of the min-plus basis terms 
used, in order to facilitate numerical computation. 

We now describe each stage in greater detail. 

A. Min-plus expansion of the terms in the value function 

This step involves obtaining a quadratic approxima- 
tion (in a min-plus sense) to the terms used in (|2]i 
(specifically (|4]i and (|5]l). The quadratic terms thus 
obtained are used to carry out the recursion step which 
yields the quadratic approximation (|6) for the value 
function at the next time step. A quadratic approximation 
of a (continuous) function g{x) over a set fi using a set 
of L quadratics is performed as follows. Note that this 
window 51 is a set centered around x* i.e., the optimal 
state estimate available at that time step. Given a set 
of points LOk G fl, we divide the region Vt into L parts 
r^i, Q.2^ ■ ■ ■ ^L- For each I e {1, 2, . . . , L} we design 
the quadratic hi{x) which approximates g{x) over Q.i 
by solving the constrained optimization problem 



there exists a set A^ -.^ J y. C and a corresponding set 
of quadratics Qm defined as 



AjJ) 



Qj + Qe, Vj e J, £eC, m(i, j) e M. 



Further the following holds 

A Qn. = A Qj + A Q^- 

D 

This result is essential to the recursion step. Specifi- 
cally, consider the following recursion equation for the 
value function (for further details/derivations c.f. [12]) 



Vi{x) 
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5(x)||dx},(9) 
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+ l\\y-c{x)\\l, (10) 



where Q is the constrained set for the choice of quadrat- 
ics defined as 



where 
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wl:=-[Q, + B^N^B]~'x[B^Lf], 
w] := -[Q„ + B'^N^B]-^ x [B'^N^Aix)], 

W := (/ + Bw\fN^{I + Bw\) + wf Q^w\. 



and subject to 
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If a discrete set of points in Vt^ are chosen in order to 
evaluate a discrete form of (|9]l, then the problem reduces 
to a least square optimization (sum of square errors). 
In specific cases there may be a simplified/efficient 
formulation to this optimization problem (this will be 
indicated for an example in Sec. llVb . 

B. Recursion to produce new quadratic approximation 
to the value function 

We start with a result on combining two sets of 
quadratics. 

Lemma 3.1: Given two sets of quadratics Qj, j E J^ 
and Qi,£ E C; For any expression of the form 



Using ([3]l, dill, Q we can write ( fTOJ i as 
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where 



A Qj + A Q^ 
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By applying Lem. 13.11 to the set of quadratics above we 
generate an index set Ii and a set of quadratics Q^^ 
such that (|6| holds. 

C. Projection/Pruning 

It can be seen that the construction of new quadratics 
leads to a growth in the number of such terms used to 
represent the value function. Hence in order to facilitate 
numerical tractability it becomes essential to reduce the 
number of such quadratics without unduly sacrificing 
estimation accuracy. This 'pruning' is, in effect, a projec- 
tion of an initially large set Q of quadratics (min-plus ba- 
sis elements) into a lower dimensional/lower cardinality 
set Qp for the min-plus expansion of the value function, 
n intuitive approach to pruning in the application of the 
max/min plus techniques in control, involves assigning a 
metric which indicates the contribution of each quadratic 
to the minimization of the value function. However in 
the filtering case such an approach can give rise to a 
non-robust filter For instance figure [Tldepicts a situation 
where the contribution of the set Q^ of quadratics to the 
minimization is much greater than that of the quadratics 
in set Qp. However in terms of robustness, if only set 
Qa was retained, it would cause the filter to respond 
more slowly to large jumps in the system state. For 
instance, an output corresponding to the state xp would 
not lead to the desired state update from Xa since all 
the quadratics around the state xp would have been 
pruned away. These situations of a jump in the state, 
arise in important applications in bimodal/multimodal 
oscillators, systems with binary or jump disturbances 
in the state. Hence an alternative pruning approach is 
required in order to increase the robustness of the filter 
One such pruning method is a clustering approach. In 
this technique we cluster (spatially) the x* the estimated 
state for each quadratic q in the original set Q. Then, 
given a requirement to retain only P of these quadratics, 
we generate P clusters and choose the quadratics with 
the best contribution metric (i.e. the ones which yield the 
most likely state) /rom within each cluster. This proceeds 
as follows 

1) for all q ^ Q, obtain the corresponding x* :— 
argminxq{x). Note; due to the presence of a 
window of interest, the argmin might differ from 
the analytically obtained minima for the quadratic 
(which might lie outside the region of interest). 

2) Now we cluster this set (say X) of the estimates 
X* such that every x £ X belongs to a cluster 
{1,2,..., P}, i.e., the we generate a cluster map 
/C which returns a cluster number for each element 
of X (based on the metric chosen to separate the 
estimates into clusters). 
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Fig. 1: Need for a clustering approach. Retaining only 
the set Qa would make it harder to identify if the state 
shifts closer to xp. 



3) Within each cluster we sort and retain the 
quadratic that has the greatest contribution to the 
minimization of the value function. Hence the set 
of quadratics to be retained is defined as 



Qp ■■= {geQ 



Vg e Q, 



/C(:e*) = lC{Xq} =^ q{Xq) < q{xq)j. 



For a visual intuition of this pruning approach, ref Fig. |2] 
This procedure yields the desired projection set Qp 
of quadratics, starting from the original set Q. 

IV. Illustrative Example 

In this section we demonstrate the ideas described 
thus far on a two dimensional system with linear dy- 
namics and a nonlinear outpuiJ. The continuous time 
representation is 



wit), 
(11) 
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■ Xl{t) ' 




'00' 




r xi{i) ] 
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dt 


_ X2(t) 
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[ ^2{t) \ 






y{t) 



40 



vit), 



where w{-) and v{-) are the process disturbance and 
measurement noise respectively. Taking a sample time 

'The implementation code for the example in this paper (also 
useable as a template for other applications) may be found at 
https://github.com/srsiidharan/robustFiltering 
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Fig. 2: The case when only 3 quadratics are to be 
retained: three clusters are formed and within each, the 
quadratic with the lowest value is selected. 



-Quadratic terms in the expansion 
-Function to fit: x^/40 
Min-plus approximation 




Fig. 3: Min-plus expansion for the output function 
a;2'^/40 for a window centered around X2 = 2. 



Af of 0.1s, the discretized dynamics is 



Xi{k 

X2{k 



1 
0.1 1 

0.1 " 




xi{k) 
■Mk) 

Aw{k), 



where Aw{k) is the approximation corresponding to 
w(fc)At over the sampling time. As specified in Sec. Hill 
the implementation of the deterministic filter proceeds 
along three steps. The general steps as modified and 
applied to this specific example are as follows to this 
example are as follows. 

1. generation of the quadratic approximation: In this 
case, the output equation ( fTTT i is 



C{x) 



x" 
40 



+ v{k). 



The contrained nonlinear optimization described in 
Sec. IIII-AI admits the following simplification in the 
current case. Consider any window J7 with subpartitions 
Q,k (recall Sec. IIII-AI ). In order to design the optimal 
quadratic approximation over such a set i7fc, we note that 
this problem is essentially a one dimensional regression 
(albeit with nonlinear constraints). We create a vector 
Xs of N sample points such that Xsi G VLk for all 
i e {1,2,...A^}. Now to fit the (continuous) output 
function (or its square), denoted by g{-), over VLk using 
a quadratic 



/(^) 
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a2 

ai/2 



ai/2 
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the optimization problem reduces to solving a least 
squares fitting problem to determine the optimal coeffi- 
cient vector z := [a2,ai, qq]'^ such that 



Az 



9{Xs2) 



where 



A 
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The additional constraints are that 02 > (for convexity) 
and that f{x) > g{x), Vx G 51. Solving this yields 
the desired quadratic bases. The min-plus basis obained 
by such a quadratic approximation over a window ft 
centered around X2 = 2 is as shown in Fig. [3] 

2. The recursion to obtain the next set of quadratics 
(i.e. the min-plus basis expansion of the value function 
during the next time step) involves taking the sum of 
the various quadratic terms in (fTOt . This is generated 
by taking the pairwise sum of the various quadratics 
available. 

3. To reduce the dimensionality of the growth in basis 
terms we use a k-means based clustering approach [14] 
where the quadratics are clustered based on the locations 
of their argmin (ref. Sec. lIII-Ct . The simulation results of 
this filter design approach for the example considered are 
as shown in Fig. IH |5] Note that although each window 
is only 2 units in width, the use of the repeated window 
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Fig. 4: State filtering 
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Fig. 5: Filtered measurement 



generation (a sliding window) leads to a filter capable 
of handling large measurement noises. 

V. Conclusion and Future Directions 

The details of the min-plus approach described herein 
generate filters for systems with nonlinear dynamics and 
nonlinear output. Its main novelty is in the utilization 
of the min-plus basis expansion of the value function 
coupled with the exploitation of the linearity of the 
dynamic programming operator over such a (semi)- 
field. The salient features of this approach as discussed 
herein, provide a structure which maybe used for a 
variety of applications. Further, the example serves as 
a template for future algorithmic developments. A few 
of the avenues along which a study of the different 
stages used in these methods may be pursued are: 
(1) generating optimal min-plus fitting techniques. The 
approach to fitting described herein provides one method 



(albeit not the optimal one); (2) the development of more 
sophisticated projection techniques for managing the 
growth in dimensionality (and an error analysis thereof); 
(3) the study of methods to speed up these stages. For 
instance some of the basis expansions may be done 
offline thereby helping speed up the real time operation 
of the implementation. The application of these methods 
to the real time estimation of signals in various domains 
is a potentially fruitful theme for future research. 
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